.. _tutorial_3:
TUTORIAL 3 : RDF
================
:Author: Andrey V. Brukhno (andrey.brukhno{at}stfc.ac.uk)
Radial Distribution Function
----------------------------
Consider a system of identical particles (atoms or molecules) at constant *N,V,T* (canonical ensemble).
The **Radial Distribution Function**, or *pair distribution function, g(r),* describes the average density variation as a function of separation between any two particles.
.. figure:: Images-RDF/RDF-illustrations-med.png
:width: 600px
To arrive at the expression for *g(r)*, we first consider the *n*-particle density distribution, which is proportional to the probablity of a *particular n*-particle configuration (:math:`{\mathbf r}_1, \dots, {\mathbf r}_n`):
.. math::
\rho^{(n)}( {\mathbf r}_1, \dots, {\mathbf r}_n ) =
\frac{N!}{(N-n)!} \int \frac{\exp[-\beta U(\mathbf r^{N})]}{Z(N,V,T)} \ d{\mathbf r}_{n+1} \dots d{\mathbf r}_N
where the *N! / (N-n)!* factor accounts for all relevant permutations (due to the particles being identical); and *Z(N,V,T)* is the canonical partition function (normalisation constant in *NVT* ensemble). Setting *n* = 2 one can define the *normalised* pair distribution function
.. math::
g({\mathbf r}_1, {\mathbf r}_2) =
\frac{ \rho^{(2)}({\mathbf r}_1, {\mathbf r}_2) }{ \rho^{(1)}({\mathbf r}_1) \rho^{(1)}({\mathbf r}_2) } =
\frac{V^2(N-1)}{N} \int \frac{\exp[-\beta U(\mathbf r^{N})]}{Z(N,V,T)} \ d{\mathbf r}_{3} \dots d{\mathbf r}_N
For an isotropic system of spherical particles one can simplify the latter formula by recognising that the *radial distribution function* depends only on the particle separation and should also be normalised by the area of the sphere of radius *r* drawn around a reference particle ("entropy on a sphere"):
.. math::
g(r)dr = \frac{V}{4\pi r^2 N^2} \langle \sum_{i\ne j} \delta(|\Delta\mathbf r_{ij}|-r) \rangle
.. g(r)dr = \frac{\rho^{(2)}(r)}{4\pi r^2 \rho_{b}} = \frac{V(N-1)}{4\pi r^2 N} \sum_{i\ne j}
.. \int \delta(|\mathbf r_i - \mathbf r_j|-r) \frac{\exp[-\beta U(\mathbf r^{N})]}{Z(N,V,T)} \ d\mathbf r^{N}
so that :math:`\rho(r)/\rho_b = g(r\to\infty)\to 1`, whereas :math:`\rho_b g(r)` gives the local density at distance :math:`r` away from a particle.
.. .. math::
.. \rho_b \int g(r) 4\pi r^2dr = N-1
NOTE: in simulation RDF can only be determined up to :math:`r_{\rm max} = \frac{1}{2} \min({\rm cell\ dimensions})`
.. figure:: Images-RDF/RDF-e1s1t12-s.png
:height: 500px
Potential of mean force (PMF)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Potential of mean force (PMF or POMF) is defined as the work done upon reversibly bringing two particles from the state of zero-interaction into the state of full-interaction at a given distance, *W(r)*.
PMF can be directly obtained from RDF:
.. math::
\beta W(r) = -\ln g(r) + \ln g(\infty)
where the zero-interaction state corresponds to infinite separation :math:`(r\to\infty)`.
NOTE: The PMF concept is also applicable to molecules, complexes or aggregates in solution.
Exercise 3.1 - RDF calculation
------------------------------
The exercises are based on :ref:`tutorial_2`, developing it somewhat further.
Change to directory **tutorial_3**, where you will find the three input files for this exercise.
Below the essential amendments in CONTROL and FIELD files for this tutorial are described. Note that you don't need to copy nor change the files from the previous tutorial!
We use reduced energy units, i.e. kT, whence Kelvin is set as energy unit in FIELD file:
.. code-block:: html
:linenos:
UNITS K
The VDW parameters are in the units of kT and Angstrom):
.. code-block:: html
:linenos:
VDW 1
LJ core LJ core lj 1.0 1.0
In CONTROL file reduced temperature is set below the critical point and pressure slightly over the critical value (note that pressure is always set in katm):
.. code-block:: html
:linenos:
temperature 1.0 # Reduced temperature T* = 1.0 < T*(CP) = 1.1876
pressure 0.02 # p* = p(katm)/0.163882576 > p*(CP)=0.1093 (reduced pressure at CP)
Collection of RDF data is invoked with the following directive in CONTROL file:
.. code-block:: html
:linenos:
sample rdf
where the three parameters are:
* **number of grid points (bins)** - set to a multiple of 10 (500)
* **cutoff** - set to half the (expected) minimum cell dimension (5.0)
* **stride (number of MC steps)** for accumulation of RDF samples - set to 1000
We also invoke collection of data for volume evolution and probability distribution, *P(V)*, with the following line:
.. code-block:: html
:linenos:
sample volume
where the two parameters are:
* **stride (number of volume steps)** for accumulation of samples - set to 10
* **bin size** - set to 10.0
Now you can re-run the *NpT* calculation for LJ system. Verify that RDFDAT.000 file has been created as a result. Note that RDFDAT file(s) are updated each time the stats are periodically saved, so one can check RDF data while the simulation is still being executed.
Plot the RDF in gnuplot::
[tutorial_3]$ gnuplot
gnuplot> plot 'RDFDAT.000' u 1:2 with lines
.. figure:: Images-RDF/Ex2-gnuplot-RDF-equil1.png
:width: 640px
Checking the results
^^^^^^^^^^^^^^^^^^^^
Questions to answer:
* Does it look reasonable?
* Any artifact observed?
You can check if the simulation converged with respect to volume (density) variations::
gnuplot> plot 'VOLDAT.000' with lines title "Volume"
.. figure:: Images-RDF/Ex2-gnuplot-Volume-equil1.png
:width: 640px
You can also check the probability distribution for volume, *P(V)*::
gnuplot> plot 'NPTDAT.000' u 1:2 with lines title "P(V)"
.. figure:: Images-RDF/Ex2-gnuplot-Vprob-equil1.png
:width: 640px
What might go wrong?
^^^^^^^^^^^^^^^^^^^^
Questions to answer:
* Has the simulation converged w.r.t. to volume?
* Are the collected statistics sufficient (noisy data, large errors)?
You might need to perform an additional run, restarting the simulation with the latest configuration stored in REVCON file.
Before doing that, store away the current results in a new directory **equil-rdf1**::
[tutorial_3]$ mkdir equil-rdf1
[tutorial_3]$ cp CON* FIELD equil-rdf1/
[tutorial_3]$ mv *.0?? equil-rdf1/
Copy the obtained RECON file into CONFIG in order to restart with the latest configuration::
[tutorial_3]$ cp equil-rdf1/REVCON.000 ./CONFIG
Based on what we learnt from the equilibration run, some amendments are necessary in CONTROL file:
.. code-block:: html
:linenos:
steps 10000000 # Number of moves to perform in simulation
equilibration 1000 # Equilibration period: statistics are gathered after this period
sample coordinate 10000 # Frequency to output configuration to trajectory file(s)
sample volume 10 2. # Sample volume: every 10th volume move (on average); bin size
sample rdf 450 4.5 1000 # Sample RDF(s): number of bins; cutoff; stride (number of steps)
After the production run and analysis, store away the results in a new directory **prod-rdf1**::
[tutorial_3]$ mkdir prod-rdf1
[tutorial_3]$ cp CON* FIELD prod-rdf1/
[tutorial_3]$ mv *.0?? prod-rdf1/
Lowering density
^^^^^^^^^^^^^^^^
In order to see the effect on RDF let us increase temperature (thereby, aiming at lower density).
Edit CONTROL file again:
.. code-block:: html
:linenos:
temperature 1.5 # Reduced temperature T* = 1.5 > T*(CP) = 1.1876
pressure 0.02 # p* = p(katm)/0.163882576 > p*(CP)=0.1093 (reduced pressure at CP)
...
sample rdf 400 8.0 1000
sample volume 10 10.
re-run the simulation and redo the analysis.
Again store away the current results in **equil-rdf2** in case you may need to redo or restart this simulation!
Consider adding the following directives in CONTROL:
.. code-block:: html
:linenos:
maxatmdist 5.0 # initial max atom displacement (set based on equilibration run)
acceptatmmoveupdate 1000 # stride (number of MC steps) for adjusting the max atom displacement
acceptatmmoveratio 0.4 # target acceptance ratio for atom moves [default: 0.37]
After the production run and analisys, store away the results in a new directory **prod-rdf2**::
[tutorial_3]$ mkdir prod-rdf2
[tutorial_3]$ cp CON* FIELD prod-rdf2/
[tutorial_3]$ mv *.0?? prod-rdf2/
Exercise 3.2 - PMF calculation
------------------------------
Calculate the PMF:s corresponding to the RDF:s obtained above::
gnuplot> plot [x=0.5:4.5] [y=-1:3] 'RDFDAT.000' u 1:(-log($2+0.000001)) t "PMF from RDF" w l
One can have both RDF and PMF in the same plot::
gnuplot> plot [x=0.5:4.5] [y=-1:3] 'RDFDAT.000' u 1:2 t "RDF for LJ" with lines, \
'' u 1:(-log($2+0.000001)) t "PMF from RDF" w l
.. figure:: Images-RDF/Ex2-gnuplot-PMF-prod2.png
:width: 640px
Further related tutorials
-------------------------
Check out :ref:`tutorial_6` and learn how to calculate free energy differences (FED).